Fit tweedie/delta models to biomass density

Author

Max Lindmark

Published

2024-02-27

Intro

Fit climate-agnostic sdms for calculating biomass trends and velocities

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "RCurl", "devtools", "modelr", "ggtext",
          "viridis", "RColorBrewer", "here", "sdmTMBextra") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)

# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# Set path
home <- here::here()

Read data

# Read data
d <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
  readr::read_csv(paste0(home, "/data/clean/catch_clean.csv")) %>%
  rename(X = x, Y = y) %>%
  pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
                 flounder_juvenile, plaice_adult, plaice_juvenile),
               names_to = "group", values_to = "density") %>% 
  separate(group, into = c("species", "life_stage"), remove = FALSE) %>% 
  drop_na(depth, temp, oxy, sal, density)
Rows: 12385 Columns: 25
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): country, haul_id, ices_rect, month_year
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)

pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, plaice_adult, …) into (group, density) [was 12385x25, now 99080x19]

drop_na: removed 9,570 rows (10%), 89,510 rows remaining
d <- d %>% filter(!species == "dab")
filter: removed 23,124 rows (26%), 66,386 rows remaining

Scale variables

d <- d %>% 
  mutate(temp_sc = as.vector(scale(temp)),
         temp_sq = temp_sc^2,
         oxy_sc = as.vector(scale(oxy)),
         oxy_sq = oxy_sc^2,
         sal_sc = as.vector(scale(sal)),
         depth_sc = as.vector(scale(depth)),
         depth_sq = depth_sc*depth_sc,
         quarter_f = as.factor(quarter),
         year_f = as.factor(year),
         year_ct = year - mean(year),
         year_sc = (year - mean(year)) / sd(year)) %>% 
  filter(quarter == 4)
mutate: new variable 'temp_sc' (double) with 12,194 unique values and 0% NA
        new variable 'temp_sq' (double) with 12,194 unique values and 0% NA
        new variable 'oxy_sc' (double) with 12,192 unique values and 0% NA
        new variable 'oxy_sq' (double) with 12,192 unique values and 0% NA
        new variable 'sal_sc' (double) with 12,195 unique values and 0% NA
        new variable 'depth_sc' (double) with 4,519 unique values and 0% NA
        new variable 'depth_sq' (double) with 4,519 unique values and 0% NA
        new variable 'quarter_f' (factor) with 2 unique values and 0% NA
        new variable 'year_f' (factor) with 29 unique values and 0% NA
        new variable 'year_ct' (double) with 29 unique values and 0% NA
        new variable 'year_sc' (double) with 29 unique values and 0% NA
filter: removed 40,246 rows (61%), 26,140 rows remaining

Read and scale the prediction grid

pred_grid <- bind_rows(readr::read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
                       readr::read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))
Rows: 457436 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 490110 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Scale variables with respect to data, not the prediction grid!
pred_grid <- pred_grid %>% 
  drop_na(depth, temp, oxy, sal) %>% 
  mutate(temp_sc = (temp - mean(d$temp)) / sd(d$temp),
         temp_sq = temp_sc^2,
         oxy_sc = (oxy - mean(d$oxy)) / sd(d$oxy),
         oxy_sq = oxy_sc^2,
         sal_sc = (sal - mean(d$sal)) / sd(d$sal),
         depth_sc = (depth - mean(d$depth)) / sd(d$depth),
         depth_sq = depth_sc*depth_sc,
         quarter_f = as.factor(quarter),
         year_ct = year - mean(d$year),
         year_sc = (year - mean(d$year)) / sd(d$year)) %>% 
  filter(quarter == 4)
drop_na: removed 4,118 rows (<1%), 943,428 rows remaining
mutate: new variable 'temp_sc' (double) with 917,493 unique values and 0% NA
        new variable 'temp_sq' (double) with 917,493 unique values and 0% NA
        new variable 'oxy_sc' (double) with 914,731 unique values and 0% NA
        new variable 'oxy_sq' (double) with 914,731 unique values and 0% NA
        new variable 'sal_sc' (double) with 897,794 unique values and 0% NA
        new variable 'depth_sc' (double) with 7,572 unique values and 0% NA
        new variable 'depth_sq' (double) with 7,572 unique values and 0% NA
        new variable 'quarter_f' (factor) with 2 unique values and 0% NA
        new variable 'year_ct' (double) with 29 unique values and 0% NA
        new variable 'year_sc' (double) with 29 unique values and 0% NA
filter: removed 471,714 rows (50%), 471,714 rows remaining

Fit the model that is preferred by most (m3, i.e., temp + temp_sq + bp_oxyg)

Use only quarter 4 since that’s what we are looking for (warmer and less oxygen than q1), and not also we use a spatial trend model Plot residuals, save model object, plot conditional effects!

# https://pbs-assess.github.io/sdmTMB/articles/spatial-trend-models.html

pred_grid_list <- list()
res_list <- list()

for(i in unique(d$group)) { 
    
  dd <- d %>%
    filter(group == i & quarter == 4)
  
  mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15)
  
  plot(mesh)
  
  # This is a variant of the main model without environmental fixed effects
  # FIXME: finalize model structure poisson-link... 
  m <- sdmTMB(density ~ 0 + year_sc + depth_sc + depth_sq,
              data = dd,
              mesh = mesh,
              family = delta_lognormal(type = "poisson-link"), # because we want to combine the spatial trends!
              spatiotemporal = "IID",
              spatial = "on",
              spatial_varying = ~0 + year_sc,
              time = "year")

  # Check the model converges
  sanity(m)
  tidy(m, "ran_pars", conf.int = TRUE, model = 1)
  #tidy(m, "ran_pars", conf.int = TRUE, model = 2)

  # Plot residuals
  # dd$res1 <- residuals(m_test3, model = 1)
  # dd$res2 <- residuals(m_test2, model = 2)
  # 
  # qqnorm(dd$res1); qqline(dd$res1)
  # qqnorm(dd$res2[is.finite(dd$res2)]); qqline(dd$res2[is.finite(dd$res2)])

  # MCMC 
  samps1 <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200, model = 1)
  samps2 <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200, model = 2)
  dd$mcmc_res1 <- as.vector(residuals(m, type = "mle-mcmc", mcmc_samples = samps1, model = 1))
  dd$mcmc_res2 <- as.vector(residuals(m, type = "mle-mcmc", mcmc_samples = samps2, model = 2))
  qqnorm(dd$mcmc_res1);qqline(dd$mcmc_res1)
  qqnorm(dd$mcmc_res2[is.finite(dd$mcmc_res2)]); qqline(dd$mcmc_res2[is.finite(dd$mcmc_res2)])
  
  # Store residuals
  res_list[[i]] <- dd
  
  # Predict on grid 
  p <- predict(m, newdata = pred_grid) %>%
    mutate(group = i)
  
  names(p)
  
  plot_map_fc +
    geom_raster(data = p, aes(X*1000, Y*1000, fill = zeta_s_year_sc1)) +
    scale_fill_gradient2()

  plot_map_fc +
    geom_raster(data = p, aes(X*1000, Y*1000, fill = zeta_s_year_sc2)) +
    scale_fill_gradient2()
    
  # Save
  pred_grid_list[[i]] <- p
    
}
filter: removed 21,377 rows (82%), 4,763 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.077508 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 775.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7571.28 seconds (Warm-up)
Chain 1:                50.626 seconds (Sampling)
Chain 1:                7621.91 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.077073 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 770.73 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7898.36 seconds (Warm-up)
Chain 1:                25.242 seconds (Sampling)
Chain 1:                7923.6 seconds (Total)
Chain 1: 

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 21,375 rows (82%), 4,765 rows remaining

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.076537 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 765.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4969.02 seconds (Warm-up)
Chain 1:                12.507 seconds (Sampling)
Chain 1:                4981.52 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.078383 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 783.83 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4781.52 seconds (Warm-up)
Chain 1:                38.625 seconds (Sampling)
Chain 1:                4820.15 seconds (Total)
Chain 1: 

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 22,374 rows (86%), 3,766 rows remaining

Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluation
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.063831 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 638.31 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4166.92 seconds (Warm-up)
Chain 1:                19.534 seconds (Sampling)
Chain 1:                4186.45 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.061829 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 618.29 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4168.19 seconds (Warm-up)
Chain 1:                19.385 seconds (Sampling)
Chain 1:                4187.58 seconds (Total)
Chain 1: 

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 22,373 rows (86%), 3,767 rows remaining

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.06067 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 606.7 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5448.28 seconds (Warm-up)
Chain 1:                20.99 seconds (Sampling)
Chain 1:                5469.27 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.064953 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 649.53 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5854.99 seconds (Warm-up)
Chain 1:                21.193 seconds (Sampling)
Chain 1:                5876.18 seconds (Total)
Chain 1: 

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 21,602 rows (83%), 4,538 rows remaining

Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluation
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.084936 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 849.36 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7757.89 seconds (Warm-up)
Chain 1:                25.68 seconds (Sampling)
Chain 1:                7783.57 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.077033 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 770.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 8032.44 seconds (Warm-up)
Chain 1:                28.351 seconds (Sampling)
Chain 1:                8060.8 seconds (Total)
Chain 1: 

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 21,599 rows (83%), 4,541 rows remaining

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.083091 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 830.91 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 9653.7 seconds (Warm-up)
Chain 1:                51.438 seconds (Sampling)
Chain 1:                9705.14 seconds (Total)
Chain 1: 
Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning: Examine the pairs() plot to diagnose sampling problems

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.079155 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 791.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 9830.03 seconds (Warm-up)
Chain 1:                52.854 seconds (Sampling)
Chain 1:                9882.89 seconds (Total)
Chain 1: 
Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: Examine the pairs() plot to diagnose sampling problems

mutate: new variable 'group' (character) with one unique value and 0% NA

# Save predictions and sims as data frames
preds_grid <- dplyr::bind_rows(pred_grid_list) %>% as.data.frame()
res <- dplyr::bind_rows(res_list) %>% as.data.frame()

Save

write_csv(res, paste0(home, "/output/residuals_04_sdms.csv"))
write_csv(preds_grid, paste0(home, "/output/data_pred_grids_04_random_sdms.csv"))

Plot residuals

# Plot residuals
# FIXME: residuals look pretty crap - delta model instead?
res %>% 
  mutate(group2 = str_replace(group, "_", " "),
         group2 = str_to_title(group2)) %>% 
  ggplot(aes(sample = mcmc_res1)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~group2) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)
mutate: new variable 'group2' (character) with 6 unique values and 0% NA

ggsave(paste0(home, "/figures/supp/qq_sdm_04_model1.pdf"), width = 17, height = 17, units = "cm")

res %>% 
  mutate(group2 = str_replace(group, "_", " "),
         group2 = str_to_title(group2)) %>% 
  ggplot(aes(sample = mcmc_res2)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~group2) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)
mutate: new variable 'group2' (character) with 6 unique values and 0% NA
Warning: Removed 12298 rows containing non-finite values (`stat_qq()`).
Warning: Removed 12298 rows containing non-finite values (`stat_qq_line()`).

ggsave(paste0(home, "/figures/supp/qq_sdm_04_model2.pdf"), width = 17, height = 17, units = "cm")
Warning: Removed 12298 rows containing non-finite values (`stat_qq()`).
Removed 12298 rows containing non-finite values (`stat_qq_line()`).